I’ve hidden all the code in this markdown by default, to make it easier to get to the map…
Continuing my explortion of Notification of Infectious Disease data, I’ve been working on an interactive map of where disease have been reported in England and Wales. I found this link a useful resource https://journocode.com/2016/01/28/your-first-choropleth-map/.
I was able to download data from the Notification of Infectious Disease website on the distrubtion of cases by local authority (“Statutory notifiable diseases - number of cases reported in week 51 of 2017”). I’ve limited this experimental analysis to one week of data. I did some refomatting in Excel and saved this as a .csv for speed and convinience, before loading it into R. I then explicitly set NA data to 0 (assuming no notifications means no cases). I also removed the string " UA" from where it appeared in local authority names.
library("reshape2")
library("magrittr")
geo_disease_data<-read.csv("geo_disease_data.csv")
geo_disease_data[is.na(geo_disease_data)] <- 0
geo_disease_data$LocalAuthority<-gsub(" UA", "", geo_disease_data$LocalAuthority)
The next task was to get some map data for local authorities. I played around with a few different maps from the Office for National Statistics Open Geography Portal (I’m new to geo-spatial data). This one seemed to be the best:
library("rgdal")
LocalAuthorities<-readOGR("https://opendata.arcgis.com/datasets/686603e943f948acaa13fb5d2b0f1275_4.geojson")
I then merged the two data sets by local authority name.
LA_Diseases<-sp::merge(LocalAuthorities, geo_disease_data, by.x="lad16nm", by.y="LocalAuthority")
The next task was to display the data as a map. For this I used leaflet. This is a wrapper for a Javascrip “App”. I’ve put some comments in the code to explain what’s going on at each step.
library("leaflet")
diseases <- names(geo_disease_data)[-1] #Get a list of disease names (the first colmn name here is local authority name)
max_cases <- max(melt(geo_disease_data, id.vars="LocalAuthority")$value) #Get the maximum number of cases for any disease
#Set up a colour scale from 0 to the maximum number of diseases
pal <- colorNumeric(
palette = "viridis",
domain = c(0,max_cases)
)
#Initialise the map with the legend and a control to allow users to select a disease
myleaflet<-leaflet(width = "100%") %>%
addProviderTiles("Esri.WorldGrayCanvas") %>%
addLegend( pal = pal,
values = 0:max_cases,
title = "Cases",
opacity = 1) %>%
addLayersControl(baseGroups=diseases,
position = "bottomleft",
options = layersControlOptions(collapsed = TRUE))
#Add a layer to the map for each disease
for (active_disease in diseases)
{
myleaflet <- myleaflet %>%
addPolygons(data=LA_Diseases,
fillColor=~pal(LA_Diseases[[active_disease]]),
fillOpacity = 0.8,
color = "black",
weight = 1,
popup = paste(LA_Diseases$lad16nm, LA_Diseases[[active_disease]],"cases"), #This popup shows the local authority name and number of cases
group = active_disease
)
}
myleaflet #show the map